% load C:\Manoj\projects\ace\20080303\season_all.mat F Txy;
% load('C:\Manoj\projects\ace\20080303\test_tf_data.mat');
%%  
%load C:\Manoj\projects\ace\20080303\test_tf_data_ap20.mat

if ispc == 1,
    load C:\Manoj\projects\ace\20080303\alldays JULI_SEG ACE_SEG N_data Txy TIME_SEG;
end;

if isunix == 1,
   load /Users/manojnair/Documents/projects/ace/20080303/alldays.mat JULI_SEG ACE_SEG N_data Txy TIME_SEG;;
   
end

%some how the first harmonic in this computation (Txy) is slightly higher
%than I do now with ace1. this creats the Tf (as created by invfrqz) with 
%tapering end (to zero)
JULI_SEG = JULI_SEG.*24.366*1e-3; %mV/m
[Txy1,F] = tfestimate(ACE_SEG,JULI_SEG,hanning(72),0,72); %Txy is agains calculated just to get
%the normalized F as required by the function invfreqz
%The TEST ACE data are corrected for 17 minutes
%[b,a]=invfreqz(Txy,F,37,0);
[b,a]=invfreqz(Txy1,F,7,4,[],3000); %better

%load C:\Manoj\projects\ace\20080303\alldays_new Txy;
%[b1,a1]=invfreqz(Txy,F,7,4,[],3000); %better

%%

if ispc == 1,
    load  c:\manoj\projects\longp\final_coefficients a_n b_n ;
end;

if isunix == 1,
    load /Users/manojnair/Documents/projects/longp/final_coefficients a_n b_n ;
end;

a=a_n;
b=b_n;


len = 120;
step = zeros([1,len]);
step(13:36) = 1; % positive  box
figure(1);
plot((1:5:len*5)./60,step,'b-','LineWidth',4);
grid on;
hold on;
figure(2);
plot((1:5:len*5)./60,filter(b,a,step).*1.0,'b-','LineWidth',4);%1.5 was for tf made only from Julia
hold on;
grid on;

step = zeros([1,len]);
step(13:end) = 1; % positive step
figure(1);
plot((1:5:len*5)./60,step,'r-','LineWidth',2);
figure(2);
plot((1:5:len*5)./60,filter(b,a,step).*1.0,'r-','LineWidth',2)

step = zeros([1,len]);
step(25) = 1; % dac delta
figure(1);
plot((1:5:len*5)./60,step,'k-','LineWidth',3);
hold on;
figure(2);
plot((1:5:len*5)./60,filter(b,a,step).*1.0,'k-','LineWidth',3)
grid on;

step = zeros([1,len]);
step(30:54) = triang(25)';
figure(1);
plot((1:5:len*5)./60,step,'c-','LineWidth',2);
hold on;
figure(2);
plot((1:5:len*5)./60,filter(b,a,step).*1.0,'c-','LineWidth',2)
grid on;


figure(1);
set(gca,'FontSize',16);
axis([0,6,-1,2]);
xlabel('Time (hours)')
ylabel('IEF Ey mV/m');


figure(2);
set(gca,'FontSize',16);
axis([0,6,-0.1,0.1]);
%set(gca,'Ytick',[-0.05, -0.025,0,0.025,0.05]);
xlabel('Time (hours)');
ylabel('Equatorial Zonal EF mV/m');

%%
% tf = conj(Txy1');
% tf2 = [tf conj(fliplr(tf(2:end-1)))];
% plot((1:5:len*5)./60, real(ifft(fft(step).*tf2)),'g-','LineWidth',2)
% grid on;



%%To compute predicted JULIA from ACE data for all selected days
load C:\Manoj\projects\ace\20080303\alldays JULI_SEG ACE_SEG N_data TIME_SEG;
load  c:\manoj\projects\longp\final_coefficients a_n b_n ;
a=a_n;
b=b_n;


for i = 1:N_data,
juli = JULI_SEG((i-1)*72+1:i*72)*24.366*1e-3;;
juli=juli-mean(juli);
ace = ACE_SEG((i-1)*72+1:i*72);
time_ax = TIME_SEG((i-1)*72+1:i*72);
pred_julia =filter(b,a,ace);
pred_julia=pred_julia-mean(pred_julia);
rms_error = sqrt(mean((pred_julia-juli).^2));
plot(time_ax+datenum(2000,1,1),juli,'r','LineWidth',2);
set(gca,'FontSize',16);
hold on;
plot(time_ax+datenum(2000,1,1),pred_julia,'b','LineWidth',2);
plot(time_ax+datenum(2000,1,1),ace./15,'k');
date_string = datestr(time_ax(1)+datenum(2000,1,1),1);
title(sprintf('date = %s, RMS err= %10.5f',date_string ,rms_error));
grid;
h=legend('Observed drift (m.s^{-1})','predicted drift (m.s^{-1})','ACE Ey mV/Km');
set(h,'FontSize',10);
%axis([-inf inf -15 15]);
datetick;
pause;
hold off;
end;

% 
%% 

load C:\Manoj\projects\ace\20080303\alldays JULI_SEG ACE_SEG N_data Txy TIME_SEG;
load c:\manoj\projects\ace\Julia_W_new w_climate Julia_W fday julia_euvac julia_solar_flux;
JULI_SEG = JULI_SEG.*24.366*1e-3;
load  c:\manoj\projects\longp\final_coefficients a_n b_n ;
a=a_n;
b=b_n;
load c:\manoj\projects\ace\OMNI_ELEC_new;
load c:\manoj\geomag\indices\aplist.mat;
ace_fday = floor(ace_all(:,1));
phase_delay = 17;
N_data=0;

for i = 1: size(w_climate,1),
    L = ace_fday == Julia_W(i).fday;
    if sum(L) > 0,
    ace_time = ace_all(L,1)+(phase_delay)/(60*24);%Try to advance 60 VS model/ ace time minutes
    ace_data = ace_all(L,2);
    L = isnan(ace_data);
    if sum(L) < 10 & sum(L) ~= 0,
        ace_data = interp1(ace_time(~L),ace_data(~L),ace_time);
    end;

    if abs(mean(diff(ace_time))-0.0035) <= 1e-004, %Use this with ACE min averages
        
        ACED = interp1(ace_time, ace_data, fday(i,:));
        pred_julia =filter(b,a,ACED)*1.0; %26-march-08 [a=7,b=4, (light tail) with filter.*1 gives relatve low rms
        %pred_julia=pred_julia-mean(pred_julia);
        pred_w(i,:) = pred_julia;
        ace_matrix(i,:) = ACED;
        missing_data_ace(i,:) = sum(isnan(ACED));
        N_data = N_data+1;
        available_ace_fday(N_data) = Julia_W(i).fday;
  
        end;
   else,
       fprintf('Day %d has some missing time stamp\n', Julia_W(i).fday);
    end;
    end;
save c:\manoj\projects\longp\pred_julia_champ pred_w available_ace_fday missing_data_ace ace_matrix;
clear;   
%% 
clear
%Calculate the RMS error between drift-model, observation etc

load c:\manoj\projects\longp\pred_julia_champ pred_w available_ace_fday missing_data_ace ;
load c:\manoj\projects\ace\Julia_W_new w w_climate Julia_W fday julia_model_drift fday;    
load c:\manoj\geomag\indices\aplist.mat;
%pred_w = pred_w.*24.366*1e-3; % pred_w is already in electric fields in
%mV/m if Txy1 is computed (in the previous cell) after scaling JULI_SEG
w_climate = w_climate.*24.366*1e-3;
w = w.*24.366*1e-3;

%size of fday is 1002, number of ace overlap is 993 - this is because ace
%does not have data on the last 9 days of julia fdays 
number_of_days = 0;


for i = 1:993,
        if length(Julia_W(i).k) >= 1,
            L = fday_ap >= fday(i,Julia_W(i).k(1)) & fday_ap <= fday(i,Julia_W(i).k(end));
            mean_ap = mean(ap(L));
        else,
            mean_ap=-1;
        end;

    if missing_data_ace(i) == 0 & available_ace_fday(i) == Julia_W(i).fday...
            & mean_ap >=20,
              number_of_days = number_of_days+1;   

         obs_climate = w(i,Julia_W(i).k);
         %obs_climate = obs_climate-nanmean(obs_climate);
         obs_minus_climate = w_climate(i,Julia_W(i).k);
         %obs_minus_climate = obs_minus_climate-nanmean(obs_minus_climate);
         obs_minus_climate_minus_ace = pred_w(i,Julia_W(i).k);
         %obs_minus_climate_minus_ace = obs_minus_climate_minus_ace-nanmean(obs_minus_climate_minus_ace);'
         climate_minus_ace = obs_minus_climate-obs_minus_climate_minus_ace;
         %climate_minus_ace = climate_minus_ace-nanmean(climate_minus_ace);

         rms_error_1(number_of_days) = (nansum(obs_minus_climate.^2));
         rms_error_2(number_of_days) = (nansum(climate_minus_ace.^2));
         rms_error_3(number_of_days) = (nansum(obs_climate.^2));
         
         var_1(number_of_days) = var(obs_minus_climate);
         var_2(number_of_days) = var(climate_minus_ace);
         var_3(number_of_days) = var(obs_climate);
         
        
    end;
end;

a=nansum(rms_error_1);
b=nansum(rms_error_2);
c=nansum(rms_error_3);
 fprintf('Number of days = %d, reduction of rms with climate = %10.5f reduction of rms climate & ace = %10.5f\n',...
     number_of_days,100*(c-a)/c,100*(c-b)/c);
 
fprintf('a = %10.5f b = %10.5f c= %10.5f\n',a,b,c);

a=nanmean(var_1);
b=nanmean(var_2);
c=nanmean(var_3);

fprintf('a = %10.5f b = %10.5f c= %10.5f\n',a,b,c);



fprintf('Ndays = %d,          model variance (observation) = %10.5f\n',number_of_days, nansum(var_3));
fprintf('             model variance (observation-climate) = %10.5f\n', nansum(var_1));
fprintf('       model variance (observation-(climate+ace)) = %10.5f\n', nansum(var_2));

fprintf('Predentage reduction climate/observations = %4.1f\n', 100*(c-a)/c);
fprintf('Predentage reduction  climate+ace/observations = %4.1f\n', 100*(c-b)/c);
%clear;
%%
%to find out the year dependency of tf
load C:\Manoj\projects\ace\20080303\alldays JULI_SEG ACE_SEG N_data Txy TIME_SEG;
for i = 1:N_data,
time_ax = TIME_SEG((i-1)*72+1:i*72);
fdayid(i) = floor(time_ax(1));
[dummy1,dummy2,dummy3] = datevec(fdayid(i)+datenum(2000,1,1));
yearid(i) =dummy1;
end;
hist(yearid,[2000,2001,2002,2003,2004,2005,2006,2007,2008]);

% this shows that number of days is equal roughly aboove and below 2003

year_limit=datenum(2003,12,31)-datenum(2000,1,1);
L = TIME_SEG <year_limit;
[Cxx,F] = mscohere(JULI_SEG(L),ACE_SEG(L),hanning(72),0,72,1/(5*60)); %1/(5*60) = sampling frequency 
figure1 = figure;
axes('Parent',figure1,'XTick',[0.1 1 10],...
    'XScale','log',...
    'XMinorTick','on');
set(gca,'FontSize',16);
box('on');
hold('all');
xlabel('period in hours');
ylabel('coherence');
hold on;
semilogx((1./(3600*F(2:end))),Cxx(2:end),'r','LineWidth',2);
axis([0.1,10,0,1]);
%found that there is no effect on cohernce 2001-2003 and 2004-2008 !



%%
for i = 1:N_data,
juli = JULI_SEG((i-1)*72+1:i*72);
juli=juli-mean(juli);
ace = ACE_SEG((i-1)*72+1:i*72);
imf = IMF_BZ_SEG((i-1)*72+1:i*72);
time_ax = TIME_SEG((i-1)*72+1:i*72);
subplot(211);
plot(time_ax+datenum(2000,1,1),juli,'r','LineWidth',2);
hold on;
plot(time_ax+datenum(2000,1,1),ace/15,'b','LineWidth',2);
set(gca,'FontSize',16);
hold on;